Noname manuscript No. 

(will be inserted by the editor) 



Momentum-Space Approach to Asymptotic Expansion for 
Stochastic Filtering and other Problems 

Masaaki Fujii 



October 4, 2012 

Abstract This paper develops an asymptotic expansion technique in momentum 
space for stochastic filtering. It is shown that Fourier transformation combined with 
a polynomial-function approximation of the nonlinear terms gives a closed recursive 
system of ordinary differential equations (ODEs) for the relevant conditional distri- 
bution. Thanks to the simplicity of the ODE system, higher order calculation can 
be performed easily. Furthermore, solving ODEs sequentially with small sub-periods 
with updated initial conditions makes it possible to implement a substepping method 
for asymptotic expansion in a numerically efficient way. This is found to improve the 
performance significantly where otherwise the approximation fails badly. The method 
is expected to provide an useful tool for more realistic financial modeling with un- 
observed parameters, and also for problems involving nonlinear measure-valued pro- 
cesses. 

Keywords Zakai equation • polynomial-function approximation ■ measure-valued 
process 

1 Introduction 

In many areas, researchers frequently encounter the situation where crucial param- 
eters for their models are not directly observable in our mother nature or in exper- 
iments. This is particularly the case, for examples, in engineering, applied physics, 
finance and economics. To get the best estimate of the unobservable from what we 
can directly observe is the goal of stochastic filtering. The most famous example with 
analytical solution is Kalman-Bucy filter [Bucy (1959), Kalman (1960), Kalman and 
Bucy (1961)], which assumes both of the signal and observation processes are linear 
and hence associated with Gaussian distribution. 

However, there are many cases where interested variables follow nonlinear stochas- 
tic processes and their distributions are far from Gaussian. This is particularly the 
case for financial problems. In fact, many people were forced to realize the sheer im- 
pacts of non-Gaussianity in the last financial crisis followed by the collapse of Lehman 
Brothers. Researchers, practitioners as well as regulators now clearly recognize the im- 
portance of understanding not only the first two moments but also every other details 
of the relevant distribution. Here, we need to deal with nonlinear filtering problems. 
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Filtering theory has a long history and is still developing very rapidly, partly helped 
by the great increase of computational power. Recently, there appeared a thick volume 
edited by Crisan and Rozovskii (2011) from Oxford university press, which contains 
latest developments and reviews for theoretical as well as numerical techniques for 
nonlinear filtering problems. Unfortunately though, they seem to require very sophis- 
ticated mathematical as well as numerical skills and there exist many hurdles to cross 
for non-specialists. 

In this paper, we propose a simple approximation scheme based on an asymp- 
totic expansion method in momentum space for nonlinear filtering problems. The 
method should be also useful for other financial problems that do not require fil- 
tering. Widely used "position-space" asymptotic expansion method (See, Takahashi 
(1999), and references therein.) is transformed into a simpler form in the momentum 
(or Fourier transformed) apace, and the resultant dynamics of the characteristic func- 
tion is given by a closed recursive system of ordinary differential equations (ODEs). It 
is shown that the form of ODEs unchanged for any order of expansion, which allows 
straightforward numerical implementation for higher order approximations. Further- 
more, dividing the original time horizon into a set of small sub-periods and solving the 
ODEs sequentially with updated initial conditions, which we call substepping method 
for asymptotic expansion, increases the parameter space where the approximation is 
effective. Two simple examples are discussed to demonstrate how the method works. 
We also make a brief comment on the possibility that the same method can be used 
to analyze other measure- valued nonlinear systems. 



2 Preliminaries for Nonlinear Filtering 

2.1 Zakai equation 

Let (f2, J 7 , P) be a probability space with a filtration (J 7 t)t>o satisfying the usual con- 
ditions. We consider n-dimensional signal process X = {X t , t > 0} and m-dimensional 
observation process Y — {Y t ,t > 0} following the dynamics of 

X t = /x(t, X t )dt + rj(t, X t )dV t + fj(t, X t )dW t (1) 
Y t = h{t,X t )dt + dW t (2) 

with Yq = and an independent initial distribution for X$. Here V and W are 
independent (P, J r )-Brownian motions with dimensionality d and m, respectively, /j,, 
h, r\ and fj are deterministic function of (t,x) and take values in R n ,PJ™,R nxd 
and R™ xm , respectively. The functions (i, r\ and fj are assumed to satisfy appropriate 
conditions so that (JTJ) has a unique solution. The measurable function h is assumed 
to satisfy the conditions that makes the following process Z = {Z t ,t > 0} be a 
martingale: 

Z t =exp^-^ h s (X s ) T dW s ~^J^ \\h s (X s )\\ 2 ds^ 

where T denotes the transposition. We denote {J^,i > 0} be the usual augmented 
filtration generated by the process Y. Our goal of the filtering problem is to obtain the 
conditional distribution 7r t of the signal X at time t given the information available 
from observing the process Y in the interval of [0,i]. In other words, for a given 
arbitrary bounded function (p, computing 



1 For simplicity, we frequently use a notation of /it(Xt) and similarly for other functions. 
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Let us define the measure P by 



Z t 



and the associated inverse relation 



Z, 



where Z+ — Z t 1 can be written as 



Z t = exp[ J h s {X s )UY s -- J \\h(X s )\\ 2 d S 



Note that the process Y becomes a standard Brownian motion in the measure P. 

We define the unnormalized conditional distribution of X to be the measure- valued 
process p = {p t , t > 0} 



Pt(<p)=E Z t <p(X t ) y t 



which is 3^-adapted and cadlag. Here, E[-] denotes the expectation in the measure . 
The desired filtered density function can then be obtained from the relation 



Mir 



It is known that the dynamics of p satisfies the following Zakai equation for arbitrary 
smooth bounded function ip: 



p t ((p) = po(<p) + / p s (A s (fi)ds+ / p s ({h s + B S ) T ipjdY t 



(3) 



with initial value po(<p) = E[(p(Xo)] associated with a given distribution of X$. Here, 
A s is the infinitesimal generator of X at time s 



n d 1 " d 2 



i dxidxj 



and 



=^(^ T (^))fc, 4 ^, k = !,■■■, m. 



For the derivation of the Zakai equation and the other technical details, see Bain and 
Crisan (2008), for example. The goal of this paper is to develop a simple scheme to 
solve the Zakai equation ([3]). 
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2.2 Filtered Characteristic Function 
Let us consider a function 

ip(£,x) = exp^ T a;) 

with £, x el", where i — V^T. If one obtains the conditional expectation of this 
function, i.e., 



for each £, it enables one to derive the conditional expectation for an arbitrary choice of 
if. This fact can be seen as follows: Let us consider the inverse Fourier transformation 

M-) 



E 


exp 




z t 


y t 


E[Z t 





M*) = 



1 



(2tt) ? 



-»£ 1 z 



which can be evaluated as 



(2tt)™ J Rd E[z t \y t ] 



E 



exp 



y 



E 


<5(*t - 






E[Z t 





The above function actually corresponds to the conditional density of the X t since 



JR d 



E 



<p{X t )Z t 



y. 



E[z t \y t ] 



Therefore, {^(^(O)} an d hence {pt(ip(£))} contains all the important information 
one needs. 



3 Asymptotic Expansion in Momentum Space 

From the discussion in the previous section, one needs to solve the Zakai equation for 
Pt(V'(0)- However, the equation looks quite complicated because of the nonlinearities 
arising from the terms p s (A s ip(£)) and p s ((h s + B s ) T i(>(£)). 



3.1 Perturbed System 

In order to make the system tractable, we introduce the perturbation parameter e and 
consider the following n-dimensional signal X^ and the m-dimensional observation 
Y processes: 

dX^ = (f t + eF t (X$ e) j)dt + (y t + e<j t (X { t e) )) dV t + e lt (X { t e) )dW t 
dY t = eH t (X ( t e) )dt + dW t 
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where now ft and v t are n-dimensional deterministic function of time and vvj is 
assumed to be positive definite 0. As before, we assume Yq = and the initial condition 
of Xq — Xq is known and given by an independent distribution. F t (x) , a t (x) , jt(x) 
and H t (x) take values in R™, R™ X<J , R nxm and R m , respectively. Furthermore, these 
functions are assumed to be given by smooth and bounded functions to guarantee the 
existence of unique solution of X^ e \ As before, V and W are independent Brownian 
motion with dimension d and m, respectively. 

As explained in the previous section, we are interested in the unnormalized distri- 
bution 



-(e)- 



with 



Z t (e) = exp (ej H s (X^) T dY s - € - J \\H s (X^)\\ 2 ds 



and 



dP 



The corresponding Zakai equation now becomes 



Here, 



3 s fe =E(%»)M^, k = l 



■ , m. 



and the infinitesimal generator is given by 

4 e) = E(^ + ei ^))^: + E 2(^ + ^W)(^ + ^( a: )) 

i= 1 * J — 1 

Our goal is to expand 

p ( t e \m) = /4 0) we)) + ep^wo) + eVfVce)) + ■ • ■ 

and obtain each p¥\?p(0) f° r J = 0> 2 • • • . 



T Q2 



ij dxidxj 



(5) 



2 We do not put e on Y since it is observed process. 
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3.2 Recursive system for Asymptotic Expansion 

We now expand the Zakai equation for each order of e. Note that, for any polynomial 
function G of x, one can write 

G(x)e* T ' = G(D i )e l ^ x 

where, G(D^) denotes the differential operator obtained by replacing each Xj in the 
function by (D^)j, which is a derivative operator defined by 

This fact allows to write 

p< e) (GWO)=GTOp<%(0) 

which is linear for p[ e \il>(£))- 

In order to avoid nonlinearity, we make use of this property of polynomial func- 
tions. With slight abuse of notations, we treat F t (x),at(x),^t(x) and H t (x) as ar- 
bitrary accurately approximated polynomial functions of x (and time) for the corre- 
sponding original functions. By Weierstrass' polynomial approximation theorem, this 
is always possible for any continuous functions within the closed interval. In practice, 
one can take wide enough interval within which the signal process resides with prob- 
ability sufficiently close to one and an associated polynomial approximation accurate 
enough for that range. 

Then, one can formally write 

4 £) m x) = (4 0) (0 + e4 X) & D t ) + e 2 4 2) (6 Dt)) x) (6) 

where 

A<?\t,D € ) = -^ T (<7 S (£> ? K T (^) + 7 ,(^)7 S T (^))C 

and similarly 

(H s (x) + B s (x)) T ^,x) = (Hj(Ds)+ieis(D e ))i>(Z,x) . (7) 

Substituting ([5]), (O and j7} into the Zakai equation (QJ, one can easily confirm the 
following result. 

Theorem 1 Each order of the asymptotic expansion p[ (ip(^)) in 

p ( t £) m)) = p { °\m) + epPm)) + e 2 P { t 2 \m) + ■■■ 

of the unnormalized filtered characteristic function 



p t (e) (^(e))=E[exp(ze T X 4 (e) ) 



00 



y, 
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satisfies 

dp[ j) m)) = A< i o \op l i j) m))dt 

+ (^ T ( J D 4 ) + ze T 7*(^))p^ 1) (^(C))^ 

wifft initial conditions p^i^iO) = Po(^(0)) Po^WKO) = ® U — 1) ^ e conven- 
tion that 

for j < 0. 

Considering a special case where there is no observation, one obtains a simple 
corollary for a standard unconditional characteristic function. 

Corollary 1 Each order of the asymptotic expansion Pt^i^iO) * n 

p { :\m) = p ( t\m) + + ^ 2) w>(0) + • • • 

of the characteristic function 

^%(£))=E[exp(^ T X t (e) )" 

satisfies 

dp t ' ) WO) = ^ 0) (OPt°' ) WO)* 



-{A«(C, J D £ )pp- 1 )(V>(0)+4 2 )(^ J D £ )p«- 2 )(V(0)}^ 



wifft initial conditions Pq^WO) = PoWKO); Po^WO) = ® (j — 1) ^ e conven- 
tion that 

p {i) m)) = o 

for j < 0. 

The above result shows that one only has to deal with a set of decoupled ODEs in 
terms of momentum {£} with a given observation path of Y. It is straightforward to 
solve the above equation for each £ up to a certain e-order, and use discrete Fourier 
transformation technique to obtain the density function. In Fourier analysis of smooth 
functions, it is well-known that most of the information is carried by small number 
of modes. In fact, in an example we provide in a later section, the resultant density 
function does not change meaningfully once the number of £-modes reaches ~ 30. This 
feature combined with the decoupled dynamics of characteristic function is expected 
to weaken the curse of dimensionality significantly, at least compared to typical PDE 
approaches. 

Analytical calculation is also possible. Since the dynamics is linear, one easily ob- 
tains the following results: 

Zeroth order 

p < i 0) (m)=eti A( ° ) ^ ds p m)) (8) 
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First order 

+ (i?7(^) +ze T 7s(^))pi 0) (^(0)^} (9) 

Higher order (j > 2) 

+ {Hj(D c ) + iei s (D s ))p^- 1 Km)dY s } . 

(10) 

Rigorous mathematical proofs of the validity and convergence of the above asymp- 
totic expansion when taking e 4- are beyond the scope of the current work. However, 
for a given observation path Y , it is likely to be proved by a similar line of argu- 
ments for the asymptotic expansion in position space without filtering problem given 
in Takahashi (1999), which is based on the results of Yoshida (1992a, 1992b, 1997) 
and Ikeda and Watanabe (1989). In the next section, we explain the inversion method 
to obtain the density function. 



3.3 Density Formula 

In this section, we provide a strategy to obtain an analytical expression of the filtered 
density. Although this is not necessary if one is only interested in numerical im- 
plementation with discrete Fourier transformation, the analytical expression is quite 
useful for various applications in finance. In particular, a model calibration and quick 
response to a client request of price indication requires very fast evaluation. 
Let us consider the inverse Fourier transformation of Pt vP {£))'• 

This corresponds to the unnormalized conditional probability density of the signal 
given the observation path of Y (See the discussion in Sec l2.2n . The desired 
normalized conditional probability density of the signal is then given by 

#(z)=4y#(*) 

where is a normalization factor defined as 

(4 e) ) _1 = f $\z)<rz. 

Thus, for applications, it suffices to calculate the expression of (z). 
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3.3.1 Gaussian distribution for Xq 

For simplicity, let us first consider the case where X is distributed by a Gaussian 
law Af(xo; Sq) with mean Xq and the covariance Sq of a symmetric positive definite 
matrix. In this case, we have 

Pom)) = I e^ Tx n[x;x ,E }d n x 



\== [ e* T *ex P (~hx-x ) T £ \x-x ))d n x 
■) n \S \ Jr" v ^ / 



where n[x; Xo, Sq] is the probability density function for a random variable with Gaus- 
sian low of M{xq; Sq), and \Sq\ denotes the determinant of Sq. The evaluation can 
be done easily by considering the variable change from x to r\ given by 

x = Xq + PqT] 

with a matrix Pq satisfying 

So = PqPo ■ 
Integration in terms of r\ easily leads to 

PoMO) = exp(^ T z - \C ZoU 
Then, from the result of previous section, we have 

(0) f.*.T 3 

r 

where 

x t = x + I f s ds 



P r=exp(it l x t --S l E t Z) (11) 



S t — S + / Vgvjds 
Jo 



Thus, it is clear that X^ has a Gaussian distribution Af(xt] S t ). If the initial position 
of Xq is exactly known as Xq = Xq, then one clearly has 

pom)) = ^ Txo 

and hence one can simply insert S = in pT|) . 

By the property of the exponential form and ©, one can check that Pt(ip(Q) is 
given by 

P?\m) = Pi°\m) (Y a s (0ds + 6.(0dy.) (12) 
with polynomial functions a s (£) € R and b s (£) € IR 1 * 1 ™ of £ 

a s (o = ^Hmr'A'HcD^p^m)) 

b s (o = p^mor 1 ( H s( D s) T +^ T is(D i )) P ^m)) 



which can be expressed by Hermite polynomials in general. 
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(27T) 

a s {D z )ds + b s (D z ) T dY s ) n[z; x t , S t ] . (13) 



t 

T, 



Thus, the first order correction to the unnormalized conditional density can be 
expressed as 



Here, 

oz 

and a s (D z ), b s (D z ) denote the derivative operator of z obtained by replacing each £ 
in the functions by D z . 

Repeating the same arguments, one can see that p^'(tp(^)) can be given, as a 



generalization of fj 12[) , by 

p[ 1 \m)=p? ) m)) f r ■■■ r(r sl ,.., Sj (Odsids 2 ---d Sj 

Jo Jo Jo v 

+r yuS2r .. >aj (£) dY Sl ds 2 ■ ■ ■ dsj 

+r yi ,y 2 ,- , yj (0 dY Sl dY S2 ■■■d), ) ■ ... (14) 

with certain polynomial functions {/ n Sli ... Sj . (£), • • • , r yu ... yVj (£)} of ^ with appropriate 
dimensions. The {• • • } denotes the term with integration order of (< j), which stems 
from the existence of the e-second order operator A^K Then, as a generalization of 
f|T3j) . the j-th order term of the unnormalized conditional density is also given by the 
correction to the Gaussian distribution: 

#W = 7Tvr / ^ t *pF(V>(0K£ 

(27T) n J R n. 

= I / •••/ (r su ... , Sj (D z )ds 1 ds 2 ■ ■ -dsj 
Jo Jo Jo v 

+r yuS2 .... , Sj (D z ) dY Sl ds 2 ■ --dsj 

+ r yi,V2,-,yj( D *) dY Sl dY S2 ---dY s ^jn[z;x t ,£ t ] H • (15) 

In the case of no (or trivial) observation, one can get the asymptotic expansion of 
unconditional probability density by putting dY terms zero. 

3.3.2 Non-Gaussian distribution for Xq 

Even when the initial distribution is not exactly Gaussian, if one can approximate it 
by the form 

4>o{z) = (some polynomial function of z) x n[z; xq, Sq] , (16) 

then the properties of the inverse transformation given in the previous section still 
hold in almost the same way. This is, for example, the case when one approximates 
the initial distribution by Gram-Charlier expansions. In the case when (|16[) holds, one 
can still write p[ ^ in the form 

p{°\i/j(£,)) — (some polynomial function of £) x exp^ T x t — — ^ T ^{^ 
and it only changes the functions {-T} in (fT4")l and (IT51) . 
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4 A Direct Application to Kushner-Stratonovich Equation 

We can also apply the technique to the Kushner-Stratonovich (KS) equation that 
describes the dynamics of the normalized conditional density of Tit instead of pt- 
Although it suffices to work on the simpler Zakai equation in filtering problems, we 
directly treat KS equation here to demonstrate the fact that the asymptotic expansion 
can also be applied to measure-valued non-linear systems. For the setup given in 
Sec l2.11 the Kushner-Stratonovich equation is given by 

ditt{<p) = TT t (A t (p)dt + {n t ({ht (<p)}(dY t -n(ht)dt) 

with a given initial value ttq (ip). This is clearly a nonlinear equation for the measure- 
valued process 7r t . See a textbook Bain and Crisan (2008) for details of the derivation. 

Let us now introduce the same perturbed system as in Sec l3.ll Then, one obtains 
the KS equation for ip(£, •) as 

+e{^ ) [(H t + B t ) T m)-4 e \H7)4 e \m)}(dY t ~e^\H t )dt) 
By the same polynomial-function approximations, one can rewrite it as 

^(#7(z? e ) + ^ T 7 *(c e ))^ . 

Now, there appears ir^ (H t ) in the equation. This term does no harm since an arbi- 
trary order j of asymptotic expansion, we need ^'(Ht) only for i = 0, 1, • • • , j — 1 
due to the additional e-factor. Thus, at the calculation of j'-th order expansion, one 
can use 

where 7Tj ("0(0) are already known for i < j — 1. 

Let us give the first few orders of expansions for the KS equation: 
Zeroth order 

d4 0) m)) = 4 a) (o4 0) m))dt (17) 

with initial value n^i'ipi!;)) = ^oi^iO)- 
First oder 

d4 1) m))=4 0) (o4 1) m))dt 

+ {(Hj(D e )+ i Cjt(D i ))n ( / , \m)-4 \Hj)^ \m)}dY t (18) 

withTT^OKOHO. 
Second order 

+ (£, D e )4 1] m)) + A{ 2) (£, (V(0) 

- (Hj(Ds) + ify t {D € ) - 4 0) (Hj))^m))4°\H t )}dt 
+{(iT7(I3 e ) + i£ T 7 t (i5 c )-7^ 

(19) 
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with initial condition 7r t (tl>(q)) = 0. 

Although the system of ODEs does not keep the same structure as the unnormal- 
ized distribution p, it is clear that one can still perform the perturbation order by 
order. Furthermore, from the discussion given in the next section, it is not always 
necessary to derive higher order asymptotic expansion for accurate estimation. 

5 Substepping Method for Asymptotic Expansion 

It is obvious by construction that the accuracy of approximation deteriorates once 
the cumulative contributions from perturbation terms 

eF, ecr, £7, eH 

become bigger. This is a common problem of asymptotic expansion methods for vari- 
ous applications. In particular for the filtering problems, requiring small perturbation 
terms seems rather restrictive since it indicates that one can treat only noisy ob- 
servations (i.e. small H). In financial applications of the position-space asymptotic 
expansion, it is known that one needs higher order approximation to reach enough 
accuracy for practical use in long-dated or high-volatility situations. There exists 
many efforts to obtain higher order corrections systematically to tackle these prob- 
lems. See, for examples, Takahashi and Takehara (2007), Takahashi et.al. (2010), and 
Li (2010) for recent developments in this direction. 

Let us consider the problem in the momentum-space approach. In Theorem [TJ we 
have seen that the relation 

dp? ) (rP(0)=Al 0) (Opi j) m))dt 

+ {A[ 1 \z,Ds)pt 1 \m)+A?\t,Ds)pt 2 \m)}dt 

+(Hj(Dt)+iejt(Ds))p i r i) (m)dY t (20) 

determines the correction terms with a given initial condition pa(ip (£,))■ Although 
higher order calculation is straightforward, there exists simpler and more efficient 
way to improve the approximation. An obvious but important feature of (|20l) is the 
fact that the recursion can be started from an arbitrary initial distribution /9o("0(£))- 
Since asymptotic expansion typically works very well for short maturities, the above 
feature naturally leads to the following substepping idea for asymptotic expansion. 

(1) Create an appropriate time grid {0 = Tq < T\ < • ■ ■ < TV = t} in such a 
way that the asymptotic expansion converges well within each sub-period [Tj_i,Tj]. 

(2) Solve h20\) of each £ for s £ [Tb,Xi] up to k-th orders of asymptotic expansions. 
This can be second (or even first) order if the stepping size is small enough. 

(3) Update the initial condition for the next period [Tj^Ta] by setting 

Pt!(^(0) *- ^eV^(^(0)j obtained in step(2) . (21) 

(4) Solve 1120)) for [Ti,T2] with the updated initial condition. 

(5) Repeat the procedures till the final period to obtain p?- vP{€)) ■ 

Although performing the above method analytically by hand is quite laborious due 



Title Suppressed Due to Excessive Length 



13 



to a large number of derivative operations, one can do it quite efficiently in numer- 
ical implementations. This is because amount of procedures required for dt and dY 
integration does not change in the above operations. In fact, one obtains accurate 
results faster by performing finer substcpping with lower-order approximation than 
performing higher-order approximation without substepping. 

The substepping method may be very useful for general measure- valued nonlinear 
equations. Although it is tedious to obtain the asymptotic expansion for complicated 
dynamics, it is definitely possible for the first few orders as we have demonstrated by 
using the Kushner-Stratonovich equation. If the approximation works well, at least 
within a very short period, the above numerical procedures will extend the effective 
region for the asymptotic expansion. If it is applied to a standard unconditional char- 
acteristic function, it should also offer an efficient option pricing method, particularly 
for long-dated and high- volatile setups. 

6 Examples 

6.1 Analytical Application to CIR Process 

Now, let us first consider the approximation of one-dimensional CIR process with no 
filtering issue, which helps to obtain a concrete image how analytical procedures work. 

We study 

dX t = 9(n - X t )dt + <Tyfx~ t dVt 

with Xq = fi. All the parameters 0, /i and a are positive constants satisfy 29 /i > 
a 1 . Then, the probability density of X t is known to have a non-central chi-squared 
distribution. 

For asymptotic expansion, we treat it as the following perturbed system: |f| 
dX ( t t] = e9F(X { t f - ] )dt + ai^ + eRiX^dVt 

with X^ = /i. Here, we have defined 
F(x) = fi — x 

R{x) = Taylor expansion at (x — /i) of (\/~x — y/JI) . 

In this example, we are going to adopt the 3rd-order expansion for R(x). Note that, 
Taylor expansion provides a good polynomial approximation only when the process 
X resides near /x. If volatility is very high, it may be better to perform a different 
method, such as the minimization of the least square differences for appropriate range. 
We shall see some examples in the next section. Systematic strategy for determining 
the optimal choice of polynomial function remains as an important future work. 
Then, the infinitesimal generator is given by 

and hence 

A^(t,Ds) = - 1 -e<J 2 R 2 m. 

From Corollary [TJ analytical calculation can be performed as follows: 

3 One needs to put e = 1 at the end for the comparison to the original model. 



11 



Masaaki Fujii 



6.1.1 Zeroth order 
We have 



with p^\ip(0) = e ' 1 ^- Thus > ii; S ives 

pi 0) Mt))= e xp(itt*-U 2 o 2 i* 



2 

Then, one obtains the zeroth order density function as 

with the definition of S t = pa 2 t. 
6.1.2 First order 

The first order correction is given by 

= f e^ 0) ^A^^D,)pf\m)ds . 
Jo 

By straightforward differential operations lead to 

= lsa 2 (8dp + a 2 )i 2 - ^isa\8p + 3sa 2 )f 
o lb 

--s pa £ + — is pa £ . 



Then, following the procedures of Sec. 13.31 one obtains 



and also 



9 . a 3 9 4 g 5 x 

' a2 d^- m3 d^ +a4 d^ +m5 d^) <t>t {z) 

with the coefficients defined by 

a 2 = -}-t 2 a 2 (89p + a 2 ), a 3 = -^-it 2 a 4 (4p + to 2 ) 
lb lb 

a 4 = -^t 3 pa 6 , a 5 = ^it^pa* . 

Higher order expressions follow similarly with the help of analytical software if nec- 
essary. In Figs. [T]and [21 we have given numerical results up to e-3rd order asymptotic 
expansions without substcpping compared with the exact non-central chi-squared 
distribution. When volatility is large, there appear sizable deviation from the correct 
distribution for small x. This is understandable because Taylor expansion near the 
origin is not accurate. Except a neighbor of the origin and (x < 0), one can see that 
our approximation well reproduces the desired density function. 
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6.2 A Numerical Application to Benes Filter 

Now, we study an example of Benes filter [Benes (1981)], where the drift of the signal 
process is nonlinear. This is a special example for which there exists an exact solution 
in the class of non- Gaussian filtering problems, and hence quite useful to test the 
current method. In the class of Benes filter, we choose a following one-dimensional 
example: 

dX t = f(X t )dt + adV t 

dY t = (htXt + h^dt + dW t 

with X = Y = 0, where f(x) is given by 

f(x) — acrtanhia— ^ . 

Here, a, <r, hi and hi are all constants. In this case, the exact filtered density of X t is 
given by 



< xact (z) 



1 , f Z\ ( h\ , w ,, A 2 



— cosh^a— J exp ( — — coth.(thi<r)z 



( r sinh(s/iicr) hi h 2 , A 

\ hl ■ uuu — \ dY s + — • r~7T7" — r COth(</licr ))z 

\ / n Sinn t/iid trsmn t/iicr a I 



Iq smh(thia) a smh(thia) a 

where n t is the normalization factor to guarantee 

7T t ° XaCt (z)dz = 1 . 

For this problem, we setup the following perturbed approximation: 

dx[ t] = eF(Xl e) )dt + adV t 
dY t = eH(xi e) )dt + dW t 

with Xg^ = Y Q = 0. Here, we use 

H(x) = h\x + hi 

F(x) — (polynomial approximation of) f(x) . 

We explain the details of polynomial approximation later. 

The infinitesimal generator contains only up to e-first order term. We have 

=eF(x)— + -a 2 — 
[) dx + 2 dx 2 



and 



A<n(£,Dt)=i£F(Dt) 



As before, one needs to put e = 1 at the end for the comparison to the original model. 
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From Theorem [1] one needs to solve the recursion for (j > 1) 

(22) 

starting from the zeroth order solution: 

Pt (0) (^(e)) = oxp(-i^v) . 

6.2.1 Polynomial- function approximation 

We now discuss how to obtain the polynomial approximation F(x) for the nonlinear 
drift of the signal: 

f(x) = atrtanh(aa;/(T) . 
Due to the normalization by a, we can roughly expect 

— ^ 1 (23) 
a 

for t £ [0,1]. This implies that Taylor expansion around x = is a natural candidate of 
F(x) when |a| < 1. When \a\ > 1, the two plateaus of f(x) start to play an important 
role in the range (|23p . Unfortunately, however, Taylor expansion does not reproduce 
the plateaus but strongly diverging behavior within the range (|23p instead, which 
destabilizes the numerical result. Thus, we take [—5a, 5a] range with a step size of 
0.2cr, and carry out least-square method (LSM) to fit a 11-dimensional odd function 
for F(x). We also adopt the weight function g{x) defined as 




with various factor w. Here, w — corresponds to a pure LSM in the 5a range and the 
polynomial function well recovers the two plateaus of f{x) in the wide range, while 
it has a relatively poor fit around the origin. On the other hand, higher w gives finer 
fit and hence finer description of the density near x — 0. In this case, however, if one 
continues to increase w it starts to destabilize the numerical result as in the case for 
Taylor expansion. Thus, we need to take a balance of this trade-off, especially when 
\a\ > 1. 

6.2.2 Numerical Results 

In the following numerical examples, we take t = 1 as the maturity and use 1, 000 
steps to create the sample observation (and signal) path. We then integrate (|2"2"j) with 
the same time step dt = 10~ 3 with a given path of Y. For differentiation, we use a 
standard finite difference method. Finally, a discrete Fourier transformation is used 
to obtain the density function. 

In the first numerical example given in Fig. [3J we have used a set of parameters 
{a = 0.8, a = 0.5, hi = 0.8, h 2 = 0.5}, and a sample path of Y given in the top graph. 
We have used w — 2.0 for getting coefficients of polynomial function F(x). The 
middle graph for the conditional density functions contains the exact one denoted by 
a red line labeled as "Benes", estimated conditional densities from (0th, 3rd, 20th)- 
order asymptotic expansion without substepping method, and those from 1st order 
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expansion with substepping method of (100, 1000) sub-periods d. One can clearly see 
the benefit of substepping method explained in Sec. [SJ Although there is no clear 
improvement from 3rd to 20th order approximation, the substepping method with 
small sub-period provides almost exact fit to the true density function. 

In Fig.m we have used hi = 10.0, which is an example of small-noise observation. 
Since a — 0.5 is relatively small, polynomial approximation for f(x) is quite accurate 
(w — 2.0). The calculation has been performed with substepping method of (100, 125, 
200, 1,000) sub-periods. Fine substepping gives almost exact density even in this case. 
In particular, the significant reduction of the variance of the conditional density due 
to high quality information provided by the observation process is well reproduced 
by the repeated application of asymptotic expansion. In this example, approximation 
without substepping is too unstable and gives only meaningless numerical results. 
As suggested in Sec. [5J the size of perturbation terms itself does not seem to be 
a relevant problem for asymptotic expansion as long as we have accurate enough 
polynomial function approximation and the substepping method. 

When one increases a, f[x) becomes like a step function and makes it difficult to 
achieve accurate polynomial approximation for the relevant range. Here, the choice of 
LSM weight w starts to affect the estimated density. In Fig.[5l we have studied the case 
of a — 1.5 with several choice of w. Here, all the calculations were done with substep- 
ping method of 1000 sub-periods. The estimated density is stable for w = 0.5 ~ 2.5, 
but becomes unstable for higher w. Note that, the impacts of LSM weight is highly 
dependent on the order of the polynomial function. As is easily guessed, the change 
of estimated density is more significant when lower order polynomials are used. In 
Fig. [51 we have done similar analysis for an example of a — 2.0, which reveals more 
clearly separated two peaks of the filtered density. 

Remark : 

In the above examples, we have used time-independent function for F(x). However, 
when there exists a significant drift for the signal process, such as large \a\ in the above 
example, making the polynomial functions time dependent is quite likely to provide 
better estimation. If we have the information about the evolution of the conditional 
mean, we can change the center of polynomial-function approximation to replicate 
the original nonlinear function more accurately in the relevant region. Initial guess 
can be obtained by extended Kalman-Bucy filter or time independent F{x) in the 
current method, for examples. 



7 Concluding Remarks 

In the paper, we have developed an asymptotic expansion technique in momen- 
tum space. Fourier transformation combined with polynomial-function approximation 
gives a closed recursive system of ODEs as an asymptotic expansion for the unnormal- 
ized conditional characteristic function. Thanks to the simplicity of the ODE system, 
higher order calculation can be performed easily. It also allows an efficient imple- 
mentation of substepping method of asymptotic expansion. As long as polynomial 
approximation of the nonlinear terms is accurate, the size of nonlinear terms ceases 
to be a big obstacle for obtaining an accurate estimation. Applications to more real- 
istic multi-dimensional filtering problems as well as other (financial) problems, such 
as option pricing (with some unobservable parameters), are left for the future research. 



5 First order expansion is good enough for short enough period since the infinitesimal generator 
contains no 2nd order term. 



18 



Masaaki Fujii 



Let us make a brief comment on the remaining problems and possible future di- 
rections of research to address these issues. As one can see, the method still suffers 
from the curse of dimensionality. However, encouragingly, there exist a large number 
of works to ameliorate the higher dimensional integration problem. See, for examples, 
Griebel and Holtz (2010), Reisinger and Wittum (2007), Reisinger and Wissmann 
(2012) and Schroder et.al. (2012). Especially, in Reisinger and Wissmann, the au- 
thors make use of the low effective dimensions of financial problems arising from a 
high correlation in the market, reduce the effective degrees of freedom significantly. 
Although they have worked in restrictive model assumptions, they succeeded to ap- 
proximate a high dimensional PDE by a series of low dimensional PDE. Applications 
and improvement of these techniques by combining the asymptotic expansion devel- 
oped in this paper looks quite interesting. 

Since we can only use finite order polynomials in practice, the quality of the 
estimation highly depends on the accuracy of polynomial approximation. When one 
has nonlinear terms difficult to fit by polynomials, the idea of change- of-variable 
developed in Takahashi and Toda (2012) may be proved to be useful. Suppose, one 
defines a new process X by using some function tf'(-) as 

X t = *{X t ) • 

If X has drift and diffusion terms that are well approximated by polynomial functions, 
one can get more accurate estimation of 

— fIp'^^ 



PM) = E e 



and hence also its density 



Then, one can recover the density of the original X t as 

where \ J\ denotes the determinant of a Jacobian matrix with the elements of (d'l r i(z)/dzj) 
Thus some of the errors can be absorbed if there exists an appropriate choice of $ '. 

Another possible solution is to use Fourier series expansion directly for nonlinear 
functions in the infinitesimal generator and observation process @. Although it effec- 
tively increases the order of integration and hence slows down the calculation, some 
functions, such as step function, are known to allow accurate approximation by rela- 
tively a small number of terms. Suppose for example, some function g has a Fourier 
expansion as 

n 

where {£ n } is a series of discretized momentum, and {g n } is a set of corresponding 
coefficients. Then, one has 



E{g(xy^}=Y.9^ e 



6 Note that, one have to resort to discrete Fourier transformation for numerical implementation 
anyway. 
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Since all the nonlinear functions are included in the perturbation terms, one can 
write the dynamics of p\ in terms of {Pt (ip(£ n )}n with k = {j — 1, j — 2}. 

As a result, it is still linear for the highest expansion order and can treat similarly as 
Theorem [T] This technique may be crucial for financial applications where stiff payoff 
functions are common. Detailed studies of these points are also among our research 
topics in the future. 
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Fig. 1 t = 3,fM = 1,<T = 0.15, 6 = 0.1. (Black, Green, Blue, Purple) lines denote (0th, 1st, 2nd, 
3rd) order approximation of asymptotic expansion, respectively. Red line denotes the exact density 
function given by a non-central chi-squared distribution. The second graph represents the difference 
from the exact density function. 




Fig. 2 t = 3, fj. = l,cr = 0.33, 9 = 0.1. (Black, Green, Blue, Purple) lines denote (0th, 1st, 2nd, 
3rd) order approximation of asymptotic expansion, respectively. Red line denotes the exact density 
function given by a non-central chi-squared distribution. The second graph represents the difference 
from the exact density function. 
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Fig. 3 t = l,dt = 10~ 3 ,a = 0.8, a = 0.5, hi = 0.8, /12 = 0.5 with polynomial function fitted with 
w = 2.0. From top to bottom, the sample path, exact and approximated density functions, and the 
difference of the approximated densities from the exact one. In the middle graph, a red line labeled 
by "Benes" denotes the exact density function. 
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— -1.20 — 

Fig. 4 t = l,dt = 10~ 3 ,a = 0.5, a = 0.5, hi = 10.0, h,2 = 0.5 with polynomial function fitted with 
ui = 2.0. From top to bottom, the sample path, exact and approximated density functions, and the 
difference of the approximated densities from the exact one. In the middle graph, a red line labeled 
by "Benes" denotes the exact density function. 
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Fig. 5 t = l,dt = 10~ 3 ,a = 1.5, a = 0.5, hi = 0.7, hi = 0.5 with 1,000 substeps with 1st order 
asymptotic expansion. From top to bottom, the sample path, exact and approximated density func- 
tions, and the difference of the approximated densities from the exact one. In the middle graph, a 
red line labeled by "Benes" denotes the exact density function. 
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Fig. 6 t = l,dt = 10 _3 ,a = 2.0,cr = 0.5, hi = 1.0,^2 = 0.5 with 1,000 substeps with 1st order 
asymptotic expansion. From top to bottom, the sample path, exact and approximated density func- 
tions, and the difference of the approximated densities from the exact one. In the middle graph, a 
red line labeled by "Benes" denotes the exact density function. 



